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We present a detailed statistical treatment of the energy calibration of hybrid air-shower detectors, 
which combine a surface detector array and a fluorescence detector, to obtain an unbiased estimate 
of the calibration curve. The special features of calibration data from air showers prevent unbiased 
results, if a standard least-squares fit is applied to the problem. We develop a general maximum- 
likelihood approach, based on the detailed statistical model, to solve the problem. Our approach 
was developed for the Pierre Auger Observatory, but the applied principles are general and can 
be transfered to other air-shower experiments, even to the cross-calibration of other observables. 
Since our general likelihood function is expensive to compute, we derive two approximations with 
significantly smaller computational cost. In the recent years both have been used to calibrate data of 
the Pierre Auger Observatory. We demonstrate that these approximations introduce negligible bias 
when they are applied to simulated toy experiments, which mimic realistic experimental conditions. 
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Figure 1. Sketch of the bias introduced by a simple energy 
cut. Shown is an ideal calibration curve S <x E of the size S of 
an air shower against its energy E (thick solid line), together 
with a low energy cut at Eq. Measured estimates of E and 
S fluctuate around the true values and spread events from an 
ideal point on the calibration curve (black dots) outwards into 
an uncertainty ellipse. The line density of events on the ellipse 
is constant. Events that migrate below the cut are discarded 
(dotted arc). Surviving events (thin solid arc) spread more 
often below the ideal calibration curve since the arc is longer. 
A least-squares fit wrongly compensates this by placing the 
calibration curve below the true curve. 


I. INTRODUCTION 

The latest generation of air-shower detectors, the 
Pierre Auger Observatory [1,2] and Telescope Array [3], 
are hybrid instruments. They combine a fluorescence 
detector which measures the calorimetric energy of an 


air shower with a low duty cycle, and a surface detec¬ 
tor array measuring its size at ground with a full duty 
cycle. 

The size of an air shower measured at the same point 
in its longitudinal development is proportional to a 
power of its energy [4], Therefore, a calibration func¬ 
tion returning an energy estimate for a measured size 
can be found by analyzing a subset of coincident events 
recorded in both detectors. 

Fitting the calibration function to pairs of energy and 
size estimates with a plain least-squares method yields 
biased results for several reasons. Firstly, the least- 
squares approach requires the true energy of the air 
shower to be known event-by-event, but the fluores¬ 
cence detector only provides an energy estimate that 
fluctuates around the true energy. Secondly, the energy 
spectrum of cosmic rays is very steep so that most of 
the data is located near the lower energy threshold of 
the detector. 

In the threshold region, the detector triggers are 
not fully efficient. Upward fluctuations have a higher 
chance of passing the trigger and entering the data set 
than downward fluctuations. This creates an accep¬ 
tance bias, so that the mean size of the accepted events 
does not reflect the true mean size of the original sam¬ 
ple. 

Applying an energy cut with a minimum energy high 
enough to avoid the threshold region altogether solves 
this problem, but it creates a new bias, caused by event 
migration over the new threshold introduced by the cut. 
How the bias appears is illustrated in Fig. 1. A superfi¬ 
cial solution is to use a slanted cut, but determining the 
angle under realistic conditions, where the resolutions 
vary with energy and size of the air shower, requires 
Monte-Carlo simulation of the data [5]. 

We will show that a probabilistic approach solves the 
problem in a consistent way. Based on the known prop¬ 
erties of air-shower development and the detectors, we 
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construct a probability density model for the experi¬ 
mental data. Maximizing the likelihood of the data un¬ 
der this model then yields an asymptotically unbiased 
estimate of the calibration curve. 


II. DEFINITION OF VARIABFES 

We use the variable S for the size of the air shower 
at the ground, where it is observed by surface detector 
arrays. The size S depends on the energy £ of the air 
shower, mass A, and geometry a. We use air-shower 
geometry as a summary term for the orientation and 
impact point of the air-shower axis. The size S is of¬ 
ten obtained by fitting an empirical lateral distribution 
function to the ground signals [6, 7], but other prox¬ 
ies work as well, such as the inferred total number of 
muons at ground in very inclined showers [8]. 

Air showers with the same geometry a and energy 
£ show a fluctuating size S at the ground. These fluc¬ 
tuations [9-11] are caused by random outcomes of the 
first few interactions of the air-shower development and 
possibly from sampling a random mass A from the 
mass-distribution of cosmic rays. The mass A is usu¬ 
ally not exactly known event-by-event and therefore the 
dependency S(A) adds to the observed fluctuations of 
S. We call these fluctuations combined intrinsic fluctua¬ 
tions. 

Our aim is to find the function that yields the mean 
size S of the air shower, averaged over intrinsic fluctua¬ 
tions, as a function of its energy £ and geometry a. The 
energy dependence is usually modeled well by a power 
law PqEP 1 . Our approach does not depend on the exact 
relationship and therefore we will just refer to p as the 
parameter vector of the function S(E, a, p). 

We mention the dependence of S on the full air- 
shower geometry a to treat the most general case. In 
practice, the dependency on a is usually corrected be¬ 
fore applying the energy calibration. The correction is 
either based on air-shower simulations [7,8], or inferred 
from data, by demanding that the flux of cosmic rays 
looks isotropic in the corrected size [12]. 

The inverse of S(E) serves as the energy calibration 
function, which provides an energy estimate £g based 
on a size S of the air shower. Care must be taken, how¬ 
ever, since the random fluctuations of the observed size 
propagate into the energy estimate. Analyses based on 
Eg need to take into account, that Eg randomly fluctu¬ 
ates around the true energy £ event-by-event, combined 
with the fact that true energies follow a very steeply 
falling distribution. This makes it more likely that a 
particular observed value of Eg was generated by an 
upward fluctuation of an air-shower of lower energy, 
than by one with the same or higher energy. If the 
distribution of energies £ is to be measured based on 
Eg [6], unfolding methods can be used [14, 15]. 

In addition to the effects discussed before, detectors 
do not measure the energy £, size S, and geometry a 


of the air shower directly. They provide estimates £, S, 
and a, that randomly fluctuate around the true values. 
These fluctuations are caused by statistical sampling of 
air-shower particles in the detector and by variations 
in the detector response. An experiment therefore pro¬ 
vides a sample of tuples (E„ S,-, a,) as input for the anal¬ 
ysis. We assume that an energy cut £ > E cut is applied 
to this set which discards events with poor resolution 
in the threshold region of the detector. 

To distinguish between functions and probability 
density functions (pdfs) in this article, we use the semi¬ 
colon in pdfs to separate the random variables from 
the dependent variables. For example, fix; p) is the 
probability density function / of the random variable 
x, whose location and shape depends on p. When inte¬ 
grals over random variables appear, we will not explic¬ 
itly indicate the limits, except if the integral does not 
cover the physical domain of the variable, for example, 
[0, oo) for £ and S. 

We will refer to the normal distributions frequently, 
and therefore use the notation AT (x; p, a) to indicate the 
density 

Ar(l; ' , '‘ ,) = 7^ exp ( _ K^) )■ (1) 

In a fully rigorous treatment, we would have to use the 
truncated normal distribution in most places, where the 
domain of the variable x is not the full real line. We 
generally assume that the experimental conditions are 
such that p/cr^> 0, so that both distributions approach 
each other. 


III. FIKEFIHOOD ESTIMATION OF THE 
CALIBRATION FUNCTION 

Our fitting method is based on the maximum- 
likelihood method [13]. For un-binned continuous data, 
it states that an estimate of the parameter vector p can 
be found by maximizing the joint pdf C of the data un¬ 
der the model considered. We make a usual substitu¬ 
tion and maximize ln£ instead of £, 

ln£(p) — ^ln/(E/, S,, a,;p), (2) 

i 

which is equivalent but easier to handle. The density 
/(E, S, a;p) models the data distribution as a function 
of p. We maximize this sum with standard numerical 
algorithms to get an estimate p of p. 

If the data density was very high, working with a 
histogram of the data would be more effective and the 
log-likelihood would take a different form. Both ap¬ 
proaches can also be combined, so that the former is 
used in high density regions to speed up the computa¬ 
tion of the sum, an example of such a technique is given 
in Ref. [15]. 
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The maximum-likelihood approach has a useful 
property that we will exploit repeatedly. Finding the 
maximum of ln£ to get the estimate p only involves the 
first derivative V p ln£. Similarly, computing the uncer¬ 
tainty estimate of p only involves the second derivative. 
Therefore, any constant factors c ; with V p C; = 0, that 
appear in the evaluation of /;(p) — /(£;, S;, ap p), can 
be dropped without changing these results, 

ln£(p) = Eln/ ; (p) = E lnc <//( P) 

i i 

= E lnc < + E ln //(p) - E ln //(p)- ( 3 ) 

i i i 

We will use this to avoid the explicit computation of 
such factors wherever possible. 

We now focus on the construction of /(£, S, a; p). The 
size function S(E, a, p) of the air-shower is at the heart 
of this pdf, the crucial point is to model the random 
fluctuations of events around this mean. 


A. Statistical model of the detection process 

The pdf /(E, S, a; p) of the observed ensemble is con¬ 
structed by folding several conditional pdfs that model 
the individual sources of fluctuations with the pdf 
h(E, a) of the arrival frequencies of air showers at the 
combined aperture of the detectors. The pdf h(E, a) 
itself is obtained by normalizing the product of the 
cosmic-ray flux /(E) = dN / (d£ dA df dQ) with the 
effective aperture A e ff(E, A, a) of the combined detec¬ 
tors. The effective aperture can dependent on energy 
and mass of the cosmic ray. For example, a fluores¬ 
cence telescope can see high-energy air showers from 
a greater distance, since they are brighter. We empha¬ 
size that for the energy calibration to work, any mass 
dependency in the effective aperture of the calibration 
data needs to be the same as in the final data set to be 
calibrated. 

We now introduce the fluctuations step-by-step. If 
the detectors were perfect and there were no intrinsic 
fluctuations, we would describe the data with distribu¬ 
tion 

/i(E, S, a;p) = S(S — S(E, a, p)) £(£, a), (4) 


placing the ^-distribution with a conditional pdf s, ob¬ 
taining 

/ 2 (£,S,a;p) = s(S;S(£,a,p),£,a) h(E, a). (5) 

The shape of the pdf s itself can depend on the air- 
shower energy £ and the geometry a. 

Now we add fluctuations caused by the detectors, 
and regard event loss from online triggers and a min¬ 
imum energy cut. These effects are modeled by an¬ 
other conditional probability density function, the de¬ 
tector kernel g(E, S, a; £, S, a). We fold / 2 with the ker¬ 
nel and multiply the result with a Fleaviside function 
0(E — £ cu t) to model the effect of the applied energy 
cut. This yields our first main result 

/ 3 (£,S,a;p) =0(£-£ cut ) J d£ J dS f da 

g{E, S, a; E, S, a) / 2 (£, S, a; p). (6) 

We neglect here, that shower-to-shower fluctuations in 
the size S may be accompanied by anti-correlated fluc¬ 
tuations in the energy estimate £. The energy estimate 
£ is typically based on light generated primarily by the 
electromagnetic cascade. A fraction of the shower-to- 
shower fluctuations is caused by variations in the flow 
of cosmic-ray energy into the hadronic and electromag¬ 
netic cascade. The anti-correlations reflect the conser¬ 
vation of energy. Anti-correlated fluctuations in the en¬ 
ergy estimate are expected to be smaller than 4 % above 
10 18 eV, and decreasing at higher energies [16]. This is 
typically small compared to the energy resolution, and 
therefore not included in the model. 

Apart from this simplification, Eq. (6) is a general sta¬ 
tistical model of the calibration data. The detector ker¬ 
nel g(E, S, a; £, S, a) can be obtained from Monte-Carlo 
simulation or derived from an empirical model. The 
pdf s needs to be estimated from air-shower simula¬ 
tions or fitted to the calibration data together with the 
function S(£, a,p). 

Due to losses modeled by the Fleaviside function and 
the detector kernel, / 3 is not normalized to unity. The 
maximum-likelihood method does not require / 3 to be 
normalized, if the normalization does not depend on 
the parameters p that are optimized. Otherwise, / 3 
needs to be replaced by 


where the Dirac ^-distribution states that the observed 
values follow the function S(E, a, p) exactly. Yet, certain 
pairs of 5 and £ occur more frequently than others, due 
to the different arrival frequencies of air showers, quan¬ 
tified by h(E, a). For a fixed geometry a, Eq. (4) repre¬ 
sents a line density embedded into the (£, S) plane. It 
traces the function S(E, a, p) that we want to extract. 

By modeling how events are randomly scattered 
away from the line density, we develop the connection 
between the function S(E, a, p) and the observed en¬ 
semble. Intrinsic fluctuations are incorporated by re¬ 


/(E,S, a;p) 


_/ 3 (£,S, a;p)_ 

JdEfdSf da/ 3 (£, S, a;p)' 


( 7 ) 


A numerical computation of the normalization is ex¬ 
pensive and should be avoided. We will show how the 
computation can be neglected in good approximation if 
a sufficiently large value of E cu t is chosen in the next 
section. 

Our approach in its general form is more complex 
than a plain least-squares fit of the data, but if it is a 
complete probabilistic model of the data, maximizing 
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the likelihood of the data under the model is guaran¬ 
teed to yield an asymptotically unbiased estimate of the 
calibration curve. In particular, the probabilistic model 
handles inefficiencies and event migration above the 
threshold defined by the energy cut £ cu t, which cannot 
be dealt with in the framework of least-squares fitting. 

The general form of Eq. (6) is expensive to compute 
due to the many integrals. We will proceed to discuss 
valid approximations which greatly reduce the compu¬ 
tational cost, up to a point of removing all integrals. 
We will illustrate these approximations along an fully 
fledged example: the application of the likelihood ap¬ 
proach to the energy calibration in the Pierre Auger Ob¬ 
servatory. 

IV. APPLICATION TO THE PIERRE AUGER 
OBSERVATORY 

Two variants of our approach have been used in re¬ 
cent analyses from the Pierre Auger Observatory. The 
variants are obtained by approximating Eq. (7) in a con¬ 
trolled manner, which reduces it into a practical form. 

Both variants have a few aspects in common. The 
Auger surface detector array is sufficiently flat and reg¬ 
ular, so that the dependence on the air-shower geome¬ 
try a reduces only to a dependence on the zenith angle 
8. The remaining atmospheric attenuation is corrected 
either empirically by demanding flux isotropy [12, 19] 
or based on air-shower simulations [18]. Thus, in this 
case the refined size parameter S depends only on the 
energy £ of the air-shower, and we have 

S(£) = p 0 £Ph (8) 

A. Detector kernel 

The detector kernel is factorized into two inde¬ 
pendent normal distributions for £ and S, and a 6- 
distribution for 8, 

g(E,S,e;E,S,e)ttg£(E-,E)g g {S-,S,e)6{0-e), (9) 

with 

y £ (E;£) =A^(£;£,ct £ (E)) (10) 

g g (S;S,8)=M(S;S,a g (S,8)). (11) 

Using independent normal distributions for £ and S 
is well motivated, since the measurements are practi¬ 
cally independent. The measurements implicitly sum 
up many individual signals with near-normal distribu¬ 
tions. Due to the central-limit theorem, the resolutions 
turn Gaussian. 

Using the ^-distribution for 0 is an approximation 
that saves a numerical integration. The ^-distribution is 
the limit of a normal distribution with vanishing width. 


so effectively this ansatz treats the measurement of 8 as 
exact. For the Pierre Auger Observatory this approx¬ 
imation is very good, since the angular resolution is 
high and the detector kernel is only a slowly varying 
function of 8. 

Event losses in the detector are neglected. This can 
be made into a good approximation by choosing a suf¬ 
ficiently high value of the energy cut E cu t. The value of 
£ cu t needs to be high enough so that all accepted events 
have 100 % detection efficiency, and rejected events with 
a reasonable chance to migrate over the threshold still 
have efficiencies near 100 %. 

Under these conditions, the normalization of Eq. (6) 
becomes independent of the choice of p. Therefore, 
Eq. (6) can be used directly in the likelihood instead of 
Eq. (7). We will drop the Heaviside function 0(E — £ cu t) 
from now on, since it is always one for the selected data 
and was only needed to compute the normalization in 
the general case. 

B. Intrinsic fluctuations 

The pdf s(S; S, £) of the intrinsic fluctuations is mod¬ 
eled by a normal distribution 

s(S;S(E,p),E,p) = Af(S;S(E,p),crs(E)). (12) 

In case of the Pierre Auger Observatory, there is no in¬ 
dication that the shower-to-shower fluctuations depend 
on the air-shower geometry a, so the dependency is 
dropped. 

The choice of a normal distribution is only empiri¬ 
cally motivated, since the pdf is theoretically unknown. 
Its shape depends on the unknown distribution of 
cosmic-ray masses. However, pure samples of simu¬ 
lated proton and iron air-showers show distributions 
close to normal [10], and the model fits data from the 
Pierre Auger Observatory well [17, 19]. 

C. Combined fluctuation model 

Since both the detector-generated fluctuations (S — 
S) and the intrinsic fluctuations (S — S) are modeled 
by normal distributions, it is tempting to carry out the 
folding 

y s (S) - J dS g g (S; S,6) s(S;S(E, p), E,6) (13) 

analytically, to save another numerical integration. Un¬ 
fortunately, the trivial solution that holds for two inde¬ 
pendent normal distributions 

g^(S) =A^(S;S(E,p),(cr? + n|) 1/2 ) (14) 

does not work here, since the resolution cr g of the detec¬ 
tor depends on the random outcome S from shower-to- 
shower fluctuations. The fluctuations are coupled. 
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which follows if the effective aperture A e g for hybrid 
events factorizes in these variables. This was found to 
hold in good approximation [10]. 

D. Variant A 

In variant A, Eq. (15) is integrated over the estimated 
zenith angle 9, which effectively removes one dimen¬ 
sion of the data from the likelihood and the inference. 
This variant has been used at the Pierre Auger Obser¬ 
vatory to calibrate air showers with zenith angles up to 
60° [19]. 

The integration leads to an effective kernel on the 
right hand side of Eq. (15) 

g?(S;E, p) - j d9g?(S;E,e,p)h e (e), (17) 


Figure 2. Exact solution g s (S;S) and normal approximation 
(S; S) to the folding of two normal distributions gg (S — S) 
and s(S — S) with crg/S = og/S = 0.2. 


The difference between the normal approximation 
gN, obtained by replacing S with S in the computation 
of erg, and the exact solution g s is illustrated in Fig. 2 for 
realistic values of crg/S — cr s /S — 0.2. The exact so¬ 
lution has a shape that resembles something between a 
normal and log-normal distribution, in agreement with 
early studies based on Monte-Carlo simulations of air 
showers [5]. 

The normal approximation is nevertheless accept¬ 
able for the estimation of S(E, p), since g^ and g s 
have the same expectation value, as pointed out in Ap¬ 
pendix A. The maximum-likelihood estimator of the 
central value of g^ is unbiased with respect to the ex¬ 
pectation value, and therefore this approximation does 
not bias the fit of the size function S(E, p). 

However, the approximation biases the estimate of 
the width <7g of the intrinsic fluctuations. While this 
parameter is irrelevant for the energy calibration itself, 
it has to be fitted to data since its true value is unknown. 
The width is also of physical interest in its own, since 
it is sensitive to the cosmic-ray mass composition. The 
impact on <7g is shown in Section V. 

By assembling all pieces together and carrying out 
the integration over 9, we arrive at the following form 

f(E,S,9; p) = J dEg i {E- l E)gf(S;E,e,p)h(E,9). (15) 

This form is the common root of two variants which 
have been independently derived and used for differ¬ 
ent zenith angle ranges at the Pierre Auger Observa¬ 
tory. Both variants assume that the distribution of hy¬ 
brid events h(E,9) factorizes, 

h(E r 9) = h E (E)he(e), (16) 


which is again approximated by a normal distribution. 
This approximation is good where the zenith angle de¬ 
pendency of the detector kernel is weak. 

The remaining partial distribution h E (E) is approxi¬ 
mated by a power law E - ' 1 ' with an appropriate spectral 
index. We arrive at variant A, 


f A (E r S;p) 



= C ( d£ —— x 
J cr E erg 

( E-E \ 2 1 / S — p 0 E pi 

V ) 2 \ erg 



E~\ 


(18) 


where frg(E) — (£rg(p 0 E Pl ) 2 + cr s (E) 2 ) 1/2 is an effec¬ 
tive resolution that includes shower-to-shower fluctu¬ 
ations, and C is an unspecified normalization constant 
that does not depend on p and is therefore irrelevant 
for the likelihood estimation. 

The log-likelihood function then is (up to constants) 


ln£ A (p) = ^ln / dE 


Cr E erg 


exp 


1 ( Ej-E 

2\ <t E 


1 fSi- 


PoE pi V 




E"T. (19) 


The integration over the true energy E is carried out 
numerically for each observed tuple (E ; -, S;). 


E. Variant B 

In variant B, the approximations after Eq. (15) take a 
different path. The zenith angle dependency is kept, 
and a bootstrap estimate is used to remove the last re¬ 
maining integral. The result is a double sum, which is 
fast to compute for small to moderate data sets. Variant 
B has been used to calibrate very inclined events with 
zenith angles beyond 65° [20]. 
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We start by observing that the pdf hf(E) of hybrid 
events is well approximated by the pdf Eg(E) in the 
range of interest, which describes the distribution of the 
observed energy estimates £. Both pdfs differ by an 
integration over the detector kernel gp of the energy 
measurement, 

hp(E) — J dEgp(E; E) h E (E). (20) 

The effect of the folding is small since the resolution of 
the energy measurement is about 10 % [2], In particular, 
in the region near E cu t and above, the pdfs differ mainly 
by a shift that can be absorbed into the normalization. 

After establishing this, we can now estimate the pdf 
h p from the calibration data itself. This could be 
done with a kernel density estimate or another non- 
parametric density estimation method. Since lip ap¬ 
pears only inside an integral, we chose the even sim¬ 
pler bootstrap estimate [21]. The bootstrap estimate 
can be regarded as a kernel density estimate with a 
kernel of vanishing width. It is constructed as a sum 
of ^-distributions positioned at the observed values £/ c , 
weighted by the inverse of the overall detection effi¬ 
ciency e k = e{E k ) at this value, 

h t^) = ^L d{E J Ek) - ( 21 ) 

K k e k 

The sum over k is independent of the sum in Eq. (2) and 
runs over all K hybrid events, not only those above the 
cut value £ cu t. 

Bootstrap theory is not well known in physics, but its 
is an established branch of statistics. For example, the 
sample mean can be derived as the bootstrap estimate 
of the expectation value for an arbitrary pdf f{x) 

E[x] = J d xxf(x) —> E b [x] = J dxx/ B (x) 

= [ dxx ^Tj S ( x ~ x k) = ^YL x k- ( 22 ) 

J K k K k 

By inserting Eq. (21) into Eq. (15) and integrating over 
E, we obtain 


/ B (E,S,0;p) = M|x 


E yM(E-,E k ,cr i/k ) hf(S-, V0 El\~a §fk ), (23) 

k tk 


with kernel width functions evaluated at the bootstrap 
values E k , 

crp k — <Jp{E k ) (24) 

1 /2 

hk = fcWf,*) 2 +1T S (E k ) 2 ) . (25) 

A final approximation replaces the output of the ker¬ 
nel width functions <7g and <Tg with their event-wise esti¬ 
mates k and iTg k from the reconstruction algorithms. 


This change is not strictly necessary, but it avoids the 
need to parametrize and fit the kernel width functions 
to data in advance. This leaves the efficiency function 
e(£) as the only external input. 

Inserting Eq. (23) into Eq. (2) yields a curious double 
sum over the energy estimates, caused by the bootstrap 
approximation. After dropping all constant factors that 
do not depend on p, which includes lig(8), we obtain 


ln£ B (p) = E ln 
i 


k ^E.k ^§,k £ k 



1 / §j ~ P0 £f 

2 V ^ S,k 



, (26) 


where the sum over i only includes events with E, > 
£cut- This is the final result of approximation B. 

The structure of Eq. (26) allows for some optimiza¬ 
tion, which may make the numerical maximization 
of the log-likelihood faster. The matrix z lk — ( E, — 
E k ) / dp k is constant with respect to changes in p, and 
can be precomputed. Terms with \zi k \ > 8 can be dis¬ 
carded altogether, as their contributions to the second 
sum are negligible. 

We note that Eq. (26) cannot be further reduced into 
an equivalent least-squares method due to the sum in¬ 
side the logarithmic terms. 


V. PERFORMANCE IN TOY SIMULATIONS 

We test the frequentistic properties of our fits on a 
set of simulated experiments. We are particularly inter¬ 
ested in the bias E[p — p] of the maximum-likelihood 
estimate p and the true coverage of the estimated con¬ 
fidence region. We approximate the confidence region 
with an ellipsoid in the usual fashion, represented by 
the covariance matrix obtained after inverting the Hesse 
matrix of — ln£(p). 

In order to study these statistical properties, we simu¬ 
late data sets of toy experiments. In these simple simu¬ 
lations data points are drawn from parametrized distri¬ 
butions, which mimic the experimental conditions for 
the detection of very inclined air showers at the Pierre 
Auger Observatory. The results obtained here hold 
equally well for a simulation of less inclined events. The 
parametrization were taken from a previous study [10] 
and are summarized in Appendix B. The distribution 
of hybrid events and the average event loss due to the 
simulated trigger is shown in Fig. 3. 

We generate 1000 toy experiments. In each toy exper¬ 
iment, we generate events until 200 pass the energy cut 
£cut = 10 18 ' 6 eV. The true size function is taken to be 


S(E) — p 0 (E/10 19 eV) Pl 


(27) 
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Figure 3. Energy distribution dN/dlnE of hybrid events 
(solid line) in arbitrary units and effective trigger efficiency 
(dashed line) used in the toy simulation. The minimum en¬ 
ergy cut used in the fits is also indicated (dotted line). 


Figure 4. Simulated data (circles with error bars) from one 
toy experiment. Also shown are the true mean size 5 (dashed 
line) of the air showers and the energy cut at 10 18 6 eV (dotted 
line). 


with input values p = (2.0,0.9). 

We fit the toy data with both variant A and B, using 
Eq. (19) and Eq. (26), respectively. For each simulated 
experiment and variant, we obtain a parameter vector 
p. After fitting many independent experiments, the av¬ 
erage (p — p) will approach the bias. 

In case of variant A, we insert the true resolution 
functions used in the toy simulation in the integrand 
of Eq. (19). The spectral index is set to 7 = —2.4, based 
on a power law fit to the simulated hybrid distribution 
between 10 18 5 eV and 10 19 ' 5 eV. 

In case of variant B, we set the efficiencies — 1 
in Eq. (26) for simplicity, since the effective trigger ef¬ 
ficiency does not depart significantly from one above 
10 18 eV. Only this region is potentially relevant for event 
migration above the energy cut. 

Intrinsic fluctuations in the toy simulation have a con¬ 
stant relative resolution, Cs/S — 0.15. To obtain results 
that generalize well, we do not assume the same in the 
likelihood fits. In principle, the relative resolution could 
vary smoothly with energy, for example, if the mass 
composition changes with energy. Therefore, we model 
the intrinsic fluctuations with a second-order Bernstein 
polynomial, 

0 's/S — qo (1 — z) 2 + £/i (1 — z)z + <72 z 2 (28) 
fO ; if £ < 10 18 eV 

z — J 1 g(E/10 llS eV) . j£ £q! 8 e y < £ < If) 20 e V 

Z 1 lg(10 20 eV/10 18 eV) ' 11 iU ev <• L Zs lu ev ' 

[l ; if £ > 10 2 °eV 

and fit the three parameters q along with the param¬ 
eters p to each data set. The parametrization with a 


Bernstein polynomial allows us to implement the phys¬ 
ical constraint cr$ > 0 with a simple lower limit qi > 0 
on the parameters, supported by most numerical opti¬ 
mization algorithms. 

While the intrinsic fluctuations are not of primary in¬ 
terest for the calibration curve, they need to be fitted in 
order to complete the probabilistic model. Assuming 
that they are zero leads to a significant bias in the fit¬ 
ted parameters p of the toy experiments. In addition, 
fitting the intrinsic fluctuations provides valuable phys¬ 
ical information, since they are sensitive to the cosmic- 
ray mass composition [ 10 , 11 ]. 

To put our likelihood methods in perspective, we also 
apply a naive least-squares fit, where p is obtained by 
minimizing 


xHp) = l( S ‘ p ° E ’’ 


S,i 


(29) 


where . is the event-wise uncertainty estimate pro¬ 
vided by the simulation. 

One of the toy experiments is shown in Fig. 4. Only 
data points above the energy cut value E cut enter the 
fits directly. In case of variant B, the data to the left is 
still used indirectly to construct the bootstrap estimate 
of /ie(E), as described in the previous section. 

A detailed comparison of the fit results for this data 
set is shown in Fig. 5, which in addition illustrates the 
dependence on the energy cut. The naive least-squares 
fit shows a large bias for any cut value in po, while 
the likelihood variants appear unbiased for energy cuts 
higher than the nominal value and give very similar 
results. A consistent bias for energy cuts below 10 18 eV 
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Figure 5. Fitted parameters for the data shown in Fig. 4 as a 
function of the energy cut value, from the likelihood fits (solid 
lines) and the naive least-squares fit (dotted line). The true 
parameter values from the simulation (dashed lines) and the 
nominal cut value are indicated (dotted vertical line). Brackets 
indicate the statistical uncertainty of the fit at the nominal cut 
value. The jitter is caused by the discreteness of the data set. 
It increases with the energy cut, since the number of accepted 
events decreases and the lever arm for the fit becomes shorter. 

_ (Po)> ((Po-po)) (Pi)’ ((Pi-Pi)) 

ML-A 2.003 ± 0.001 (+0.003) 0.898 ± 0.001 (-0.002) 

ML-B 1.998 ± 0.001 (-0.002) 0.898 ± 0.001 (-0.002) 

Naive LSQ 1.910 + 0.001 (-0.090) 0.892 + 0.001 (-0.008) 

Table I. Averages of the fitted parameters from 1000 toy exper¬ 
iments for the two likelihood fits and the naive least-squares 
fit. The estimated bias is shown in parentheses. 


is also found for the likelihood fits. The jitter in the 
scans varies from experiment to experiment, but these 
biases appear consistently. 

The bias in case of low energy cuts is expected, since 
the trigger efficiency begins to deviates significantly 
from one below 10 18 eV, as shown in Fig. 3, and the se¬ 
lection bias due to the trigger becomes important. Since 
efficiency terms were taken out of the statistical model 
to obtain the approximations, the fits are not applicable 
at such low energy cuts. 

The biases are further explored in Fig. 6 and Table I. 
The naive least-squares fit is strongly biased, while the 
observed bias in the two likelihood fits is negligible. 
The statistical uncertainties in a single toy experiment 
are an order of magnitude larger. The estimated confi¬ 
dence regions constructed from log-likelihood cover the 
true values in 66 % of the toy experiments for variant A 
and B, very close to the expected 68 %. 

The fitted intrinsic fluctuations <7g are shown in Fig. 7. 


Figure 6. Averaged parameter estimates from 1000 toy ex¬ 
periments for the two likelihood fits (circles) and the naive 
least-squares fit (diamond), compared to the true values (dot¬ 
ted lines). Shown in the background are the individual out¬ 
comes obtained from variant B (small dots), those of variant A 
are very similar. The inset zooms closer to the true values. El¬ 
lipses around the average fitted values represent the statistical 
uncertainty of the average. 



lg(E/eV) 

Figure 7. Averaged estimate of the intrinsic fluctuations eg / S 
from 1000 toy experiments for the two likelihood fits (thick 
lines with error bands), compared to the true constant value 
(horizonal dots). Shown in the background are the individual 
outcomes obtained with variant B (thin lines), those of variant 
A are very similar. 
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Both likelihood fits show a small bias /S — cr s /S) of 
about 10 % below 10 19 eV. At higher energies, variant 
B performs slightly better than variant A. The varia¬ 
tion in the fit from experiment to experiment is large, 
however, especially at the high end of the energy range. 
The scatter at the high end reflects that the fit is less 
constrained where the data density is low. The inter¬ 
pretation of trends in the fitted fluctuations therefore 
has to be done carefully. 

Further exploration showed that the bias vanishes, if 
the toy simulation is done with the normal approxima¬ 
tion to the combined fluctuation model instead of 
the correct computation. The bias is therefore a conse¬ 
quence of the normal approximation described in Sec¬ 
tion IV C. 

In a final test, we push the bootstrap estimate of 
/*g(E) in variant B to the extreme by generating another 
set of 1000 toy experiments with only 10 events above 
E C ut- Again, we find only negligible bias in the param¬ 
eters of the calibration curve. 


VI. CONCLUSION 

We presented a statistical model that describes data 
taken by a hybrid air-shower detector, consisting of a 
surface detector array, which measures a size of an air 
shower, and a fluorescence detector, which measures an 
energy estimator. We developed a maximum-likelihood 
approach based on the statistical model, which allows 
us to infer an asymptotically unbiased energy calibra¬ 


tion curve for the size from coincident events observed 
in both detector parts. 

Since the general model is somewhat cumbersome to 
use, we derived two approximations. The approxima¬ 
tions lead to handy formulas, without sacrificing ac¬ 
curacy or introducing significant bias. Both approxi¬ 
mations are used by the Pierre Auger Observatory in 
different zenith angle ranges. 

We applied the more aggressive approximation to 
simulated toy experiments, to investigate the statisti¬ 
cal performance. The results showed that the estimated 
calibration curves are not biased with respect to the true 
curve used in the generation of the toy data. 
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Appendix A: Expectation of ^(S;S) andg s (S;S) 

We compute the expectation value of Eq. (13) and 
Eq. (14). The experimental conditions that we consider 
are §/ u[S] 0 and S/ tr[S] 0. This allows us to com¬ 
pute the expectation value over the whole domain of 
real numbers in good approximation. The expectation 
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value for Eq. (14) immediately follows 

/ 0° 

dSSg*(S;S) = S, (Al) 

-oo 

since y^(S;S) is a normal distribution around S. 

For the computation of the expectation value of 
Eq. (13), we use that by definition the expectation val¬ 
ues of gg(S;S) and s(S;S) are S and S, respectively By 
changing the order of integration, we obtain 


The reconstruction codes in the real experiment pro¬ 
vide event-wise estimates of these resolutions, which 
are used by fits. The estimates randomly vary around 
the true resolutions. To simulate this, we multiply the 
true resolutions that are used internally with a factor 
(1 + O.lz) before they are passed on to the fits, where z 
is standard normal-distributed. 

The reduced trigger efficiency for small air-shower 
sizes is included in the toy simulation. The trigger 
probability has the form of a cumulative normal dis¬ 
tribution. 


E[S] = / dS S / dSg s -(S;S)s(S;S) 
J — OO J — 00 

/ oo _ pOO 

dSs(S;S) / dSSy g (S;S) 

-OO J—oo 

poo 

= / dSs(S; S) S — S. 


(A2) 

P(S, 9) = 

(A3) 

2 ^(0) J 


}i{6) = (l-z)po + zpi, z = 

(A4) 

cr(0) = (1 -z)p 2 +Zp 3 , Z = 


p = (-0.95,-1.3,0.2,0.6). 


0/0.35, 

0/0.35, 


Appendix B: Toy simulation of calibration data 

We make a Monte-Carlo simulation of Eq. ( 6 ) to ob¬ 
tain artificial calibration data. We start by simulat¬ 
ing the arrival distribution h(E,9) — He(E) hg(9) of air 
showers at the detector aperture. The events generated 
here simulate highly-inclined events with 60° < 0 < 
80°, as they are seen by the Pierre Auger Observatory. 
In the numerical formulas that follow, the energy £ is 
in units of eV, the zenith angle 0 in units of radian. 

The energy spectrum li£ is modeled by a broken 
power law with a low-energy suppression that takes the 
form of the cumulative normal distribution. The zenith- 
angle spectrum hg is modeled as modified exponential 
decay. 


The probability is taken to be a function of the size es¬ 
timate S, not S, since the trigger decision is highly cor¬ 
related with the sampling fluctuations of S around the 
true size S. The reduced trigger efficiency of the en¬ 
ergy measurement is effectively included in the energy 
distribution /?e(E) and not simulated separately. 


hE(E,c ' eik ( j£ ^r) x 

( EP2-0.3 ; if 17.0 <lg£ < 18.3, 
x < E ? 2 ; if 18.3 < lg E < 19.6, 

[ePz- 12 ; if 19.6 < lgE, 

hg(9) cx exp(p 3 z + P 4 Z 2 ), z = 0 - 1.047, 
p = (18.3,0.3,-2.3,-6.4, -45.0). 

The shower-to-shower fluctuations are taken to be 
normal-distributed and with a constant relative reso¬ 
lution cr s /S — 0.15. 

The detector-generated fluctuations for the energy es¬ 
timate £ and the size estimate S are drawn from normal 
distributions with relative resolutions, 

a i /E=l P0 + Pl(1§£ “ 18 ' 4)2 ; if 1 § £ ^ 18 ' 4 ' 

1 po ; otherwise, 

<T S /S — p 2 + P3/VS, 

p = (0.10,0.03,0.04,0.10). 




